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Abstract 

We solve a physically significant extension of a classic problem in the theory of diffusion, 
namely the Ornstein-Uhlenbeck process [G. E. Ornstein and L. S. Uhlenbeck, Phys. Rev. 
36, 823, (1930)]. Our generalised Ornstein-Uhlenbeck systems include a force which de- 
pends upon the position of the particle, as well as upon time. They exhibit anomalous 
diffusion at short times, and non-Maxwellian velocity distributions in equilibrium. Two 
approaches are used. Some statistics are obtained from a closed-form expression for the 
propagator of the Fokker-Planck equation for the case where the particle is initially at rest. 
In the general case we use spectral decomposition of a Fokker-Planck equation, employing 
nonlinear creation and annihilation operators to generate the spectrum which consists of 
two staggered ladders. 
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1 Introduction 



This paper introduces a physically important extension of a classic problem in the theory 
of diffusion, namely the Ornstein-Uhlenbeck process jjj. Our results are obtained by 
spectral decomposition of a linear operator. The spectrum of this operator consists of two 
ladders of eigenvalues with respectively odd and even parity. The ladders of eigenvalues 
are staggered, that is the odd-even step is different from the even-odd step (see Fig. 
The corresponding eigenfunctions are generated by a raising operator. A concise account 
of our work on these staggered ladder spectra appeared earlier j^j- In the following we 
show how the results summarised in |2j were obtained. We also derive new results, not 
included in our earlier report: a closed-form solution for example, and the generalisation 
of our previous results to a continuous family of diffusion processes. 

1.1 The Ornstein-Uhlenbeck process 

Before we discuss our extension of the Ornstein-Uhlenbeck process, we describe its usual 
form 1 . This considers a particle of momentum p subjected to a rapidly fluctuating 
random force f(t) and subject to a drag force — 7P, so that the equation of motion is 

P = -7P + /(*)■ (1) 

The random force has statistics (f(t)) = 0, (f(t)f(t')) = C(t — t') (angular brackets 
denote ensemble averages throughout). If the correlation time r of f(t) is sufficiently 
short (7T <C 1), the equation of motion may be approximated by a Langevin equation: 

dp = — 7£>di + dw , (2) 

where the Brownian increment dw has statistics (dw) = and (dw 2 ) = 2Dodt. The 
diffusion constant is 

Do = \ j°°Jt(f(t)f(Q)). (3) 
This problem is discussed in many textbooks (for example [2]). 

1.2 Generalised Ornstein-Uhlenbeck processes 

Our extension arises when the force depends upon position as well as time. We consider 
the case where the fluctuations of the force on the particle are mainly a consequence of 
the spatial, rather the temporal, fluctuations of the force f(x,t). A consequence of this 
difference is that the impulse dw supplied to the particle in a short time St depends upon 
the momentum of the particle. If the particle is at position xq at time to, this impulse is 

rt +6t 

5w = dtf(x +p(t-t )/m,t)+O(St 2 ). (4) 

Jto 

In particular, the impulse approaches zero as the speed \p\/m of the particle increases, 
because the motion of the particle effects an average over the spatial fluctuations of the 
force. This can be seen clearly by considering the second moment of Sw. We assume that 
the force f(x,t) has the following statistics 

(f(x, t))=0, (f(x, t)f(x', t')) = C(x -x',t- t') . (5) 
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The spatial and temporal correlation scales of the random force f(x,t) are £ and r respec- 
tively. We consider the case where (for most of the time) the momentum of the particle is 
large compared to po = m£/r, then the force experienced by the particle decorrelates more 
rapidly than the force experienced by a stationary particle. If 5t is large compared to r 
but small compared to I/7, we can estimate the variance of the impulse (5w 2 ) = 2D{p)5t 
as follows (due to translational invariance, we consider without loss of generality a particle 
which starts from position x = at time t = 0) 

{5w 2 ) = / dh dt 2 (f(pti/m,ti)f(pt2/m,t 2 )) 
Jo Jo 

/oo 
dt C(pt/m,t) + 0(T 2 ) . (6) 
-00 

We define the momentum diffusion constant by writing (5w 2 ) = 2D(p)5t + 0(5t 2 ), and 
find 

1 f 00 

D{p) = - dtC(pt/m,t). (7) 
^ J -00 

When p <C po we recover D(p) = Dq. When p 3> po, we can approximate to obtain 

Dip) = + 0(p- 2 ) , D X = ^ [°°dX C(X, 0) . (8) 
\P\ 2po J -00 

When the force is the gradient of a potential V(x, t) with a correlation function having 
continuous derivatives, we find that D\ is zero. This case is discussed in section |SJ where 
it is shown that D(p) ~ \p\~ 3 provided the correlation function of V(x,t) is sufficiently 
differentiable. Another variation, also discussed in section |HJ arises when the correlation 
function of the force exhibits a discontinuity at t = (as when the potential V(x, t) is 
itself generated by an Ornstein-Uhlenbeck process). In this case D(p) ~ \p\~ 2 , and other 
exponents are also possible. We therefore consider a general situation where D(p) ~ 
and give exact results for case 

D(p) = D c (p Q /\p\)( (9) 

with £ > 0. We analyse the dynamics by solving a Fokker-Planck equation which de- 
termines the probability density for the Langevin process in which the momentum has 
diffusion constant given by ©• We discuss the form of this Fokker-Planck equation in 
section [5J the remainder of this introduction will set our work in context with earlier 
research on related topics. 



1.3 Earlier work 

The motion of a damped particle subjected to a force fluctuating in both space and time 
was first studied by Deutsch [1], who addressed an entirely different aspect of the problem. 
Deutsch considered to case where the momentum of the particle remains small compared 
to po, and posed the question of whether particles aggregate. He discovered that there is a 
phase transition between coalescing and non-coalescing trajectories. (Two of the authors 
of the present paper subsequently solved Deutsch's one-dimensional model exactly 0, and 
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results for two and three spatial dimensions are discussed in |H1 El ) • All of these papers 
only considered cases where p <C pq. 

Sturrock [S] analysed the motion of a particle subjected to a spatially varying force field 
without damping. He introduced the concept of a momentum diffusion constant which 
varies as a function of the momentum: that is, he considered the same problem as is 
addressed in the present paper, but in the limit of damping constant 7 = 0. Subsequently 
Golubovic, Feng, and Zeng |5j identified the importance of the relation D(p) ~ \p\~ 3 (in 
the case of a potential force), and discussed the nature of the Fokker-Planck equation and 
its solution in the case where 7 = 0. It was argued that the particle exhibits anomalous 
diffusion and solution for the propagator of the Fokker-Planck equation with initial value 
p = was proposed. Later Rosenbluth ^0] pointed to an error in the evaluation of this 
propagator. The results of [SI El DDI were applied to the stochastic acceleration of 
particles in plasmas, and subsequent contributions have concentrated on refining models 
for the calculation of D(p) (see, for example, |121 113j ). 

In the following we analyse the problem with the damping term, proportional to 7, 
included. Surprisingly, we find that this more general problem is more tractable: we are 
able, for example, to obtain precise results concerning the problems considered in 
by taking, in our solutions, the limit 7 — > 0. 

There is a large literature devoted to the motion of particles advected in random 
velocity fields (corresponding to the large-7 limit of the model we study). In the case where 
the velocity field is independent of time, sub-diffusive motion is typically found [TH IT5] . 
The advection of tracers in a turbulent fluid is described by models with rapidly fluctuating 
velocity fields |16j . In our problem the inertia of the particles plays an important role. 
Particles suspended in a turbulent fluid can show surprising clustering properties when 
inertia effects are significant. These were first proposed by Maxey |17j : the current state 
of knowledge is summarised in 0. In cases where the random force results from motion 
of the surrounding fluid, it is not possible for the condition p S> po to be realised [7j- 

We remark that a brief summary of many of the results of this paper was already 
published 2 . The closed- form solution of section |3J the WKB analysis and most of the 
results for general values of £ were not discussed in this earlier work. 

1.4 Description of our results and outline of this paper 

In order to simplify the presentation, we describe in detail only our results for the case 
C = 1, corresponding to a generic random force. Corresponding expressions for general 
values of C are obtained using the same method, and we quote the most important results 
for general values of £ in section |S] at the end of the paper. 

In section |2 the Fokker-Planck equation for the generalised Ornstein-Uhlenbeck pro- 
cesses is described. In section 01 we briefly discuss a particular closed-form solution, 
which enables us to determine the steady state momentum distribution (which is non- 
Maxwellian) and some statistics, such as the time evolution of the variance of the mo- 
mentum. The results of section El are not sufficient to enable all statistics to be calcu- 
lated, and in the general case we obtain statistics via a spectral decomposition of the 
Fokker-Planck equation. Section 0] discusses this spectral decomposition. We transform 
the Fokker-Planck operator into a Hermitean operator and determine the eigenvalues and 
eigenvectors of this 'Hamiltonian' operator by generating them using a new type of raising 
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and lowering operators, which are nonlinear second-order differential operators. We show 
that the resulting spectrum is a ladder spectrum, consisting of separate ladders for the 
odd and even parity states. These are staggered: the odd-even separation differs from 
even-odd. Section |S] contains calculations of the matrix elements needed for computing 
correlation functions and expectation values. In section we summarise our results on 
diffusion and anomalous diffusion for generic random forcing. 

Section [7| discusses a technical issue concerning our evaluation of the spectrum. When 
the index of the eigenvalue is large, it is possible to apply standard WKB approximation 
methods everywhere except in the vicinity of a singularity of the Hamiltonian. We show 
that the singularity introduces phase shifts which explain the staggered-ladder structure 
of the spectrum. 

Finally, in section |H1 we explain in more detail how other values of £ can arise and 
summarise our results for general £. 



2 Fokker- Planck equations 

We consider a particle with equations of motion 

x=p/m, p = — r fp + f(x,t) (10) 

where the force f(x,t) is random, with statistics given by equation In the limit as the 
correlation time r of the force approaches zero, the equation of motion of the momentum 
may be approximated by a Langevin equation, (p|). where the random increment dw has 
second moment (dw 2 ) = 2D(p)dt with D(p) given by ([7jl. This Langevin equation for the 
stochastic evolution of p(t) corresponds to a Fokker-Planck equation (generalised diffusion 
equation) for the probability density of the momentum, P(p,t). Using standard results 
[3], the Fokker-Planck equation is 

f = -|(»w^) + ^(^) ("J 

Note that we can replace dp by dw in the expression for D(p), because the neglected terms 
are of higher order in dt, and that D(p) has already been obtained in equation (JJJ). In 
order to determine the correct form of the Fokker-Planck equation it remains to determine 
(dp) = - jdt + (dw). 

Expanding the impulse (@J) about a reference trajectory x(t) = pt/m, and using the 
fact that (f(x,t)) = 0, we obtain 

(5w) = ^Jdttjdhjtoa exp[-7(fe-«3)]^drfi/m,ti)/(pi 3 /m,*3)V(13) 

Note that throughout the three-dimensional region of integration, we have < £3 < 
t-2 < *i < St, and the short correlation time implies that the integrand is negligible 
unless \t\ — t$\ < t. The integrand is therefore significant along a line rather than a 
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surface, because t% must lie between t\ and £3. The integral is therefore 0(5t), rather than 
0(5t 2 ) which would obtain if the integrand were significant on a surface. We replace the 
factor exp[— j(t 2 — £3)] by unity because 7T <C 1, and the other factor is negligible when 
1*2 — ^3 1 > t. The integral over ti then gives simply t\ — t^. Writing t = ti — £3, in the 
limit jt <C 1, the result is therefore 

{sw) = ^SZ dtt {^ Mf{pt,m,t) ) =5t ^p D{p) - (14) 

This implies 

vip) = ^D(p) . (15) 

Rosenbluth jlUj has pointed out that this relation can also be obtained as a consequence 
of applying the principle of detailed balance. 

With ((7J) and (|15|) . the following Fokker-Planck equation obtains: 

Sturrock 8 introduced a related Fokker-Planck equation (without the damping term) and 
also gave an expression for D(p) analogous to equation 0. 

In the following we discuss our solution of (|16|) for the particular case of generic random 
forcing (corresponding to £ = 1). Results for other values of £ are obtained in an analogous 
fashion. The general case is briefly described in section |SJ 



3 A particular closed-form solution 

In this section we introduce a particular solution of the Fokker-Planck equation (|16|) 
with D(p) given by (jHJ). We restrict ourselves to the case of generic random forcing 
(corresponding to £ = 1): 



— = —( +D ——)P 

dt dp \p\ dp 



(17) 



Consider the distribution P(p,t) of momentum p for particles initially at rest. It satisfies 
the initial condition P(p, 0) = 5(p) where 5(p) is the ^-function. For this particular initial 
condition, we have found the following closed-form solution of (|17j) : 



,1/3 



2r(4/3) [3p £>i(l - e- 3 -r*)]V3 



exp 



l\p\ 



3p £>i(l - e" 



(18) 



Eq. ()18|) determines how the moments of momentum grow for a particle initially at rest: 
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r(i/3) 



(19) 



for positive integers /. This result is consistent with the result obtained in [5] (eq. (8) in 
that paper). In the limit of small times (|19|) gives rise to anomalous diffusion 



(20) 
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At large times {yt S> 1), by contrast, we obtain a stationary non-Maxwellian momentum 
distribution 

1 T 1 / 3 

PM ~ 2T(4/3) (3p,g, )'/» ■» I " *I 3 /< 3 P° D '>] ' < 21 > 

The particular solution 1)18(1 generalises in a natural way to other values of £. 

However, in order to determine the momentum correlation function and the spatial 
diffusion properties, the particular solution (|18jl is not sufficient, the general solution for 
arbitrary initial condition is required. We have not been able to obtain the general solu- 
tion to 1)16)1 in closed form. Therefore, we determine it using spectral decomposition: we 
construct the eigenvalues A n and eigenfunctions ip n (z) of a Hermitian operator H corre- 
sponding to the Fokker-Planck equation l|17j). We identify raising and lowering operators 
A + and A which map one eigenfunction to another with respectively two more or two 
fewer nodes. We use these to obtain the spectrum and eigenfunctions of H which in turn 
allow us to construct the propagator, expectation values and correlation functions. This 
approach is described in sections to |SJ 



4 Spectral decomposition 

Introducing dimensionless variables (t 1 = jt and p = zpo[D% / (^Pq)] 1 ^ 3 ) we write ()17)l as 

BP 8 Id 

% = w^ + v\i )PsFP (22> 

It is convenient to transform the Fokker-Planck operator Ftoa Hermitian form which we 
shall refer to as the Hamiltonian operator: 

H = P - 1/2 FP 1/2 = - - ^ + ^ A# • (23) 
2 4 dz\z\dz v ; 

Here Pq(z) oc exp(— |z| 3 /3) is the stationary solution 1)21)1 satisfying FPq = 0. We solve the 
diffusion problem by constructing the eigenfunctions of the Hamiltonian operator. In the 
following we make use of Dirac notation |18j of quantum mechanics to write the equations 
in a compact form and to emphasise their structure. 

The eigenfunctions of the Fokker-Planck equation 1)161) are alternately even and odd 
functions, defined on the interval (— oo,oo). The operator H, describing the limiting case 
of this Fokker-Planck operator, is singular at z = 0. We identify two eigenfunctions of 
H by inspection, iPq(z) = Cq exp(— |z| 3 /6) which has eigenvalue Aq = and ipQ (z) = 
Cq z\z \ exp(— |z| 3 /6) with eigenvalue Aq = —2. These eigenfunctions are of even and odd 
parity, respectively (zero and one node, respectively). Our approach to determining the 
full spectrum is to define a raising operator A + which maps any eigenfunction (z) to 
its successor with the same parity, ^ +1 {z), having two additional nodes. 

4.1 Algebra of raising and lowering operators 

We write 

H = a~\z\- l a + . (24) 
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Here = (d z ± z\z\/2). We introduce the operators 

A = a + \z\~ 1 a + and A+ = ar\z[~ x ar (25) 

as well as 

G = a + \z\- 1 a-. (26) 
Note that A + is the Hermitian conjugate of A. The commutator of A and A + is 

[A,A+} = -3(H + G). (27) 

Note also that H — G = I (where / is the identity operator). 

4.1.1 Eigenvalues 

It can be verified that 

[Ft, A] = 3A and [H, A + ] = -3i + . (28) 

These expressions show that the action of A and A + on any eigenfunction is to produce 
another eigenfunction with eigenvalue increased or decreased by three, or else to produce 
a function which is identically zero. The operator A + adds two nodes, and repeated action 
of A + on iPq(z) and i^q{z) therefore exhausts the set of eigenfunctions. Together with 
= and Aq = —2 this establishes that the spectrum of H is (see Fig. ^) 

= — 3n and A~ = — 3n — 2 n = 0, . . . , oo . (29) 

4.1.2 Eigenfunctions 

We represent the eigenfunctions by of H by kets |^~) and ). The actions of A and A + 
are 

A + \i>t) = C± +1 |^ +1 ) and %±) = C± +1 h/,±_i>- (30) 
The normalisation factor C~ +1 is determined as follows: 

1 = (VVhIC+i) = (C- +1 )- 2 ^-\AA + \^) = (C- +I r 2 (i;-\[A,A + ] +A+A\^) . (31) 
It follows 

(Cn+i) 2 = (3 (-2A- + 1) + {C-f) . (32) 

By recursion we obtain 

C~ = [3n(3n + 2)] 1 / 2 . (33) 
This determines the normalisation of the states (A + ) n \ipQ) 

\i>-) = N-(A+) n \^) (34) 

with 

(n \ - 1 / 2 

3fe(3fc + 2) J ]V -. (35) 



For the positive-parity states we proceed in a similar fashion and obtain 

C+ = [3n(3n-2)] 1 / 2 . (36) 

This implies 

\ri)=N+(A+r\rt) (37) 

with 

(n \ - 1 / 2 

n3fc(3fc-2)J iV+. (38) 

The operators A + and A differ from the usual examples of raising and lowering operators 
in that they are of second order in d/dz, whereas other examples of raising and lowering 
operators are of first order in the derivative. The difference is associated with the fact that 
the spectrum is a staggered ladder: only states of the same parity have equal spacing, so 
that the raising and lowering operators must preserve the odd-even parity. This suggests 
replacing a first-order operator which increases the quantum number (total number of 
nodes) by one with a second-order operator which increases the quantum number by two, 
preserving parity. 

There is an alternative approach to generating the eigenfunctions of H. This equation 
falls into one of the classes considered in [19] . and we have written down first-order oper- 
ators which map one eigenfunction into another. However, these operators are themselves 
functions of the quantum number n, making the algebra cumbersome. The approach is 
briefly described in the next section. 



4.2 Schrodinger factorisation 

Consider the eigenvalue problem 

m) = m (39) 

with Hamiltonian (|23j) . For z > it can be transformed by the variable change x = z 3 : 



A 2 

dx 2 



_d_ 

dx 



(3x) 2 — r + 3x- x[ —x + A — - 







(40) 



by means of the transformation x = z 3 . Eq. (|4()j) is a Fuchsian linear differential equation 
with regular singular points of rank less than or equal to two. Eq. Q4U|) can therefore be 
factorised using a generalised Schrodinger factorisation scheme (see |19| for a review of 
this method). 

Applying this scheme we have obtained raising and lowering operators generating the 
spectrum (|29|l. The raising operator acting on ip^ is given by 



T ±in+1 = -3x-^ + | - 3(n + 5±) 
with 5+ = 1/3 and 6- = 1. The lowering operator acting on ip^ is 

T±, n = -3x-^- - \x + 3(n + 77±) 
dx 2 



(41) 



(42) 
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with ?7 + = and t?_ = 2/3. Note that the Hermitian conjugates of 2± n +i and T-j- n are 



(T ± , n+1 ) + = -T ±i „ +1 + 3 (43) 

(f±, n )+ = -T±, n + 3. (44) 

The raising and lowering operators satisfy 

T ±!n+1 |^±) = C± +1 |V± +1 ) and f ± , n |^) = C±|^±_ 1 ) (45) 



with given by (j33j) and (|36|) . generating the spectrum (|29j) . The operators differ from 
A and A + introduced in section^in that they are of first order in d/dx, and in that they 
depend on the state they are applied to. 

5 Correlation functions and matrix elements 
5.1 Correlation functions 

The required solutions of the Fokker-Planck equation may be expressed in terms of the 
propagator K(y, z, t) which is the probability density for the scaled momentum to reach z 



after time t, starting from y. It satisfies the Fokker-Planck equation dfK = FK and can 

— 1/2 

be expressed in terms of the eigenvalues A^ and eigenfunctions <ftn(z) = P Q 4 , n( z ) °f F'- 

K(y, z; t') = <(yWn(z) exp(A^) (46) 

na 



for t' > 0. Here y is the initial value and z is the final value of the coordinate. The 
expansion coefficients a^(y) are determined by the initial condition K(y, z\ 0) = 5(z — y), 
namely a°(y) = Pq ij}^(y). In terms of the eigenfunctions of H we have (for t > 0) 

K(y,z;t>) = ^P - 1/2 (y)C(y)P 1/2 (z)CWexp(A^) . (47) 

■na 

Equilibrium correlation functions of an observable O(z) are given by 

/oo poo 
dz dyO(z)O(y)K(y,z;t')P (y). (48) 
-oo J — oo 

Since P^ 2 (y) = ipgiy) this corresponds to 

(O(z )O(z t ,)) cq . =J2\(^\0\r n )\ 2 eMKt') (49) 

na 

for t' > 0. The momentum correlation function in equilibrium, for instance, is 

/ n \ 2/3 

(pt P0 ) eq . = Pi ^ Y.^t\^n? exp(A-t') (50) 
\1PqJ 

for t' > 0, which requires the evaluation of matrix elements (iPq\z\i/j^). 
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Consider on the other hand the time-dependence of (x 2 (t)), with particles initially at 
rest at the origin. We need to evaluate 



77!/ 



t ft 



(x 2 {t)) = \ I dti / dt 2 ( Ptl pt 2 ). (51) 



Jo 



In dimensionless variables this corresponds to 

<* 2 (t)> = i (to*) 2/3 -L r r df 2 <^> . (52) 

7 V 7 / m Jo Jo 12 
The required correlation function is (assuming i 2 > ij. > 0) 

/OO /'OO 
dzi / dz 2 ziz 2 K(zi,z 2 ; t' 2 - t[) K(0, Zi;t[) 
-oo J — oo 

= E^^o + l^l^)(^I^IV'+>exp[A-(t 2 -t' 1 ) + A+t' 1 ]. (53) 

In order to evaluate (|53|). the ratios of wave-function amplitudes ^+(0)/^>g"(0) are required 
in addition to matrix elements of z. The matrix elements (V'n l-^IV'rn) are determined in 
section 15.21 while the ratios of eigenfunctions are calculated in section 15.31 

5.2 Matrix elements 

To evaluate the matrix elements Z mn = {ip^\z\tp^) we proceed in three steps: we first 
evaluate Zq u , in a second step, the matrix elements Z mn are related to Zo in _ m for m < n. 
Third, Z mn is evaluated for m > n. 

5.2.1 Matrix elements (tp^ \z\ip^) 

Consider first Zq u = (tPq\z\i/j~). These matrix elements are obtained by recursion. To 
evaluate 

Z , n+ i = (i>+\zA + \4,-)/C- +l (54) 
we write zA + = zG + z{A+ - G) = z{H — J) + z(A+ - G). It follows 

(V + |M + |^-) = (A- - l)Z 0n + (^\z(A + - G)\^~) . (55) 
Using (A + — G) = —za~ and [ii 2 ,a~] = —2z we obtain 

(4;+\zA + \4;-) = (\; i +l)Z 0n . (56) 
This corresponds to the recursion 

zon = {-it , n^Qfc + i) Zqo (57) 

With Z 00 = 3-^ 12 y/¥/T(2/3) (found by direct evaluation of an integral) we obtain 
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5.2.2 Matrix elements (ip^ n \z\-tp n ) for m <n 

Consider now the case m < n. Let 

Jmn = W+\A m z(A+r\^) = <^ + |i m M + ](i + ) n -Vo~> + (^\A m A+z{A+) n - 1 \^) ■ 

We use [z, A + \ = — (|z| _1 a _ + d - !^" 1 ) to write 

Jmn = -(V'o + ii m Nr 1 a-(i + ) n ~ 1 iV'o~)-(V'o + i^ m ^i^r 1 (^ + ) n " 1 iV'o") 

+(^\A m A + z(A + r- 1 \^) 

= /(I) + 7(2) , 7(3) (KQ) 
"mn ' "mn 1 " mn \'-" J ) 

and evaluate the three terms separately. The third one gives 

jffi = (^ + |i m i + i(i + )"-Vo ) = (C^f^A^KA+Y-^) = (C+) 2 J m _ ln _x . 

(60) 

Consider next the first term: A m \z\~ x ar (A + ) n - 1 = i m - 1 a+|z|" 1 d+|z|" 1 d-(i+) n - 1 = 
A m - 1 a + \z\- 1 G(A + ) n - 1 , and thus 

^+|l m |i|- 1 o-(l+)^ 1 |^o> = (K-i- 1 )(^\A^ 1 & + \z\-HA+)^ 1 \^). (6i) 

Using A m - 1 a + \z\- 1 {A + ) n - 1 = A^a^z^a-^^cT {A + ) n ~ 2 = A m ~ l G\z\- 1 aT A +n ~ 2 
we obtain the recursion 

4& = (^_i-l)(A+_ 1 -l)ji 1 l lfl _ 1 = 3n(3m-2)J r S- ln _ 1 (62) 
Now consider J^. Using a^\z\^ l (A + ) n ^ 1 = A + \z\~ l a~ A +n ~ 2 it follows 

Jmn = (^m) ^m-ln-1 • (63) 

This implies 

7( 2 ) = ( C ^) 2 7(1) = ™ /(I) (RA\ 

u mn iV\ + i s mn "mn V"^V 

Note also that J® = 0, as well as Jq„ = 0. This gives Jq* = J^ n [consistent with (|64|)]. 
We obtain 

m 

^n!n = II 3 ( n - m + A: )( 3fc - 2 ) Jo «- m ( 65 ) 
fc=l 

which results in 

7(D - ( ^m+n r(2/3)3 m +"+ 5 / 6 r(n + l)r(m + l/3)r(n - m + 1/3) 
Jm " ~ 1 lj 6^ r(n-m + l) ' [ ' 

Eq. (|6U|) allows us to write down an inhomogeneous recursion for J mn : 

Jmn = 3m(3m - 2) J m -ln-l + (l + m/n)J« (67) 

12 



where the inhomogeneous term is given by (|66|) . Iterating this recursion we obtain 



m I m 



J mn = E( II j (l+_i_)jg)^ H (68) 

Z=0 \fc=2+l / 

which results in 

_ w 3 m +"+ 1 r(2/3) 2 r(n + l)r(m + l/3)r(n-m+l/3) 

J TO „ = (-1) —2 m+n+1 — — — J 00 . (69) 

47r z 1 [n — m + 2) 

This determines J mn for n > m. 

5.2.3 Matrix elements (V'ml^lV'ri) f° r m > n 

For m > n we use instead 

= <^|i m ii +n IV> ~) = (Vtfl^^^+IV'o"} + (^ + l^ m " 1 ^(^ + ) n I^X70) 

and proceed as before. We find that J mn = for m > n + 1. For m = n + 1 we obtain 

Jn+m = {^W\- l a + {A+r\^) + {^\A n a + \z\-\A 

- jM +J(2) +J(3) f711 

~~ ^n+ln ~ J n+ln ~ "'n+ln V ' A J 

Consider first J^+ ln : 

Jn+ln = (Cfi ) Jnn-l = 3n(3rt + 2)J n J_ 1 . (72) 



Second, using A n |z| 1 d + (A + ) n = A n 1 a + |z| 1 J 4( J 4 + ) ra we determine J^lin' 

■ff+m = (Cn ) 2 ^li = 3n(3n + 2)^ . (73) 

Third, 

Jnlm = (A+-l)(A-_ 1 -l)^l 1 = (3n + l)3n (74) 

We deduce that 



_ 3n + 2 f(2) 

and obtain the recursion 



Jn+ln = J { nl ln + J^ ln + J^ ln = + 2)J nn _ x + (l + . (76) 

Iterating (|74j) we obtain 



■Sin = (II 3*(3fc + l)) ^ = ^9^(2/3^ + 1^ + 4/3)700 (77) 



v,fc=l 
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and thus from (|76|) 



/3 

Jn+m = ^-9 n r(2/3)r(2 + n)T(n + 4/3) J 00 • (78) 



Comparing this result to (|fi9j) we find that gives the correct result for m = n + 1, 
although it was derived assuming m < n. Normalising to obtain Z mn our final result is 

Zmn = Mr-fV + „ + 1 )r(2/3) V£^^^ r '"-"' + '/f (T9) 

6vr v / r(m+ l)r(n + 5/3) T(n - m + 2) 

for n > m — 1 and zero otherwise. 

5.3 Ratios of eigenfunctions 

In this section we show how to evaluate tp^ (0)/V ; o"(0)- For z > 0, the eigenfunctions Vvt" (2) 
are of the form N^g n (z) exp(— z 3 /6), where g n (z) is polynomial in z 3 , of the form 

9n{z)=g^ +g^z 3 + ... (80) 

We determine how A + and J? act on these polynomials. To this end we define 

i'+ = ^/^A + e- zS ' & = {d z -z 2 )z- 1 {d z -z 2 ) (81) 
= e z3 / 6 He~ z3 / 6 = (d z -z 2 )z- 1 d z . (82) 

This implies 

i'+ - = -(<9 Z - 2 2 )z = z 3 - zd z - 1 (83) 

and thus 

{X + -H')g n = -g^+0{z 3 ). (84) 

Using = — 3n we obtain 

= -3ngW + 0(z 3 ) (85) 
from the eigenvalue equation. Taking (|84|) and l|85|) together we have 

i'+ 5n = _(3 n + l) 5 (0) +C , (z 3). (86 ) 

Eq. (|5H|) implies 

^ = -(3^ + 1)^°). (87) 

With (1301) it follows 

V++i(0) = iV+ +1 & = -(3n+l^ (88) 



Our final result for the ratio of wave-function amplitudes is therefore 

Vvt(o)/V + (o) = (-i) 



v / 3r(2/3)r(n + l/3) 



2vr T(n + 1) 
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6 Equilibrium correlations, diffusion and anomalous diffu- 
sion 



In this section the momentum correlation function in equilibrium and the time dependence 
of {x 2 (t)) are determined. 

6.1 Momentum correlation function in equilibrium 

The correlation function of momentum in equilibrium is obtained from (|49jl. We have 

(PtPo}eq. = Pi — Yl Z 0n eX P( A n *') ■ ( 90 ) 

V7Po/ n 

Using (|79*|l we obtain 

for t > 0. Here F 2 i is a hypergeometric function [2U]. It follows that (ptPo)eq. decays as 
exp(— 2~/t) at large times as opposed to exp(— jt) in the Ornstein-Uhlenbeck process. 

6.2 Diffusion at long times 

We now turn to (x 2 (t)}. This expectation value is calculated using eqs. (|52|). (|53|) . (|79j) . 
and (JHU). We have 

<x> W ) _ («^) 2,3 4 T f f (92) 



with 



Jo Jt' x Jo Jo 

-e^)-Ar(l- 



A+(l-e A r*') _A-(l-e A ^ i ') , . 

2— • (93) 



We define 



a = ft^ 7 7 _ 3 2 / 3 r(2/3) 2 (fc+/+i)r(fc+i/3)r(/+i/3)r(/-fc+i/3) 

H ~^(0) ' 12vr 2 T(k + l)T(l + 5/3)T(l-k + 2) { > 

for > I — 1 and zero otherwise. In order to determine (x 2 (t)}, the sum 

S{t') = Y J MiT k i{t') (95) 

kl 

is required. Note that A^i = for k < I — 1. Consider the behaviour of (|93|) at large values 
of i'. We write the k = term separately 

T« = — + -^(e^'-l). (96) 

\. \ 
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This gives 

(x 2 (t)) = 2V x t (97) 

/Po-Di \ 2 / 3 1 A 2 2 Kt , 

+ I 7 J mV^ 0! A- 2( j 



+ (— — ") o 2 ^ ^ TTTT^ Z 0l Z klTkl(t') . 



At large i' the secular term dominates and diffusion is thus recovered. The diffusion 
constant is obtained as |2j 

_ (poDx) 2 / 3 ^3~ 5 / 6 dlH 5 x 
x m 2 7 5/3 2r(2/3) 2 32 V3'3'3'3'3' /• 1 j 

6.3 Anomalous diffusion at short times 

Next we consider short times. In order to evaluate (|95|) at small values of t', we replace 
the sums in (|93|) by integrals: 

roo pi 

S(t')= dl dk T(k,l,t)A(k,l) . (99) 



J o jo 

In order to evaluate ()99() we use the asymptotic form of the coefficients A^: 

3 2 /3r(2/3) 2 k + j 

127T 2 ^!4/3(!-^/3' ^ 

The coefficients A(/c, I) exhibit non-integrable divergence k —* I. In view of this divergence 
we make use of a sum rule of the A^: 

£Ajw = 0. (101) 

fe=o 

It can be derived by considering 

!-f 1 oo oo 

^twwtiwr) = 5>* (*)^k wr> = E<*W>^fc = <*wr> (i° 2 ) 

fc=0 fc=0 fc=0 

which vanishes for z = 0. Replacing sums by integrals the sum rule amounts to 

f dk A(k,l) = 0. (103) 
J o 

Eq. ()1U3() allows us to write 

coo /■( 



/'OO W 

S(t') = / dl dfc 
io Jo 



T(k,l;t')-limT(k,l;t') 
k—>i 
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A(Jfe,l)- (104) 



We find that the divergence of A(k, I) is reduced to an integrable divergence by the fact 
that T(k, I; t) — limfc_>; T(k, I; t) = 0{k — I). Approximately, T(k, I; t) is given by 



T(k,l;t) 



3fc[l - exp(-3Z*)] - 3Z[1 - exp(-3fct)] 
27lk(l - k) 



(105) 



Changing the integration variables in (|1U4|) to x = 3lt and xy = Skt, we have k = xy/(3t) 
and I = x/(3t), and 



127T 2 2/2/3(1 -y)5/3' 

The Jacobian of the transformation is J = x/(9i 2 ). In the new variables, 



(106) 



T(x,y;t) 



-ail 



a(x) — a(xy) 



where a(x) = [1 — exp(— x)]/x. Furthermore 



lim T(x, y; t) = —2 — lim 



a(x) — a(xy) 

i^ 



Using 



we finally obtain 



lim T(x, y; i) = —2 — 
S/-+1 a; 



d_ 

dy 



a(xy) 



y=l 



-2— xa'(x) 
x 



T(k,l;t) - limT(k,l;t) 

k— >l 



This gives 

poo f'l 

S(t) = -Ct 8/3 / dx x- 8/3 / dy 
Jo Jo 



a(x) — a(xy) 



xa [xj 



a(x) — a(xy) 



xa (x) 



1 + 2/ 



y 2 /3(l - y)5/3 



(107) 



(108) 



(109) 



(110) 



(111) 



with C = 3 1 / 3 r(2/3) 2 /(27r 2 ). Equation (|111|) implies anomalous spatial diffusion at short 
times: 

<x 2 (t)) = C x {(poD^m- 2 ) t 8 / 3 (112) 

with 



f'OO f'l 

C x = -C / dx x~ 8 / 3 / dy 
Jo Jo 



a(x) — a(xy) 



xa'(x) 



1 + y 



y 2/3(! -2/)5/3 ' 



(113) 



This anomalous diffusion is analogous to that described by Golubovic, Feng and Zeng 
in the case where there is no damping (they considered the case where £ = 3 in their 
paper). Our results give the prefactor as well as the scaling behaviour. It is noteworthy 
that we are able to solve the problem discussed in |5] by solving a more complex set of 
equations exactly, and taking a limit. 
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7 WKB analysis 



The staggered ladder spectrum discussed in section|I]is surprising, especially in view of the 
fact that for large quantum number n we expect that the eigenfunctions of the Hamiltonian 
(|23|) might be obtained by WKB theory. In this section we show how to obtain the 
spectrum (j29j) of H using asymptotic Wentzel-Kramers-Brillouin (WKB) analysis. We 
show how phase shifts associated with the singularity at z = of the Hamiltonian (|23|) 
are the source of the staggered spectrum. It turns out that the WKB procedure gives rise 
to the exact eigenvalues. 

In dimensionless coordinates, the classical Hamilton function corresponding to (|23j) is 



H d = l -- Z ^-p 2 /z. (114) 

Solving H c \ = A for p we obtain p(z, A) = ±|\/ (2 — 4X)z — z 4 , while the velocity is 

z = dH cl /dp = -2p/z . (115) 

The classical trajectories are figure-of-eight orbits, illustrated in figure^ 
The WKB wavefunction is of the form 

f(z) = [z/p(z, A)] 1/2 exp Ui J\\zp{z, A)) . (116) 

The phase of the wave function can be determined as follows. We discuss separately 
the behaviours of the wavefunction at the origin zq = and in the vicinity of the regular 
turning point zt. p . = (2 — 4A) 1 / 3 . In the latter case, the wavefunction at z < z t . p . connected 
with the turning point is 

ft. P .(z) = C t . p .[z/p(z,X)} 1/2 sm (J* tP 'dzp(z,\) + ^ . (117) 

Consider now the behaviour of the wavefunction near the origin. The Hamiltonian 
has a singularity at z = 0. We find an exact solution of the equation when 7 = 0. This 
equation has a continuous spectrum, but we identify solutions f+{z) and f-(z) which 
correspond respectively to even and odd solutions of the full equation for 7 > (with 
discrete spectrum). Close to z = 0, the damping term in the Hamiltonian is negligible, 
and for z > the eigenfunctions resemble solutions of the equation 

d z z- l 8 z f{z) = -Af(z) (118) 

where A(= —A) is a positive constant. Write / = F' , and find that 

so that F" + AzF = Cz for some constant C. Thus we find G(z) = F{z) — C/A satisfies 
G" + zAG = 0, which has solution G(z) = Ai(— A 1 / 3 ^), and a similar solution constructed 
from Bi(x) (here Ai(y) and Bi(y) are Airy Ai and Bi functions [201 )• ^ ne general solution 
is 

f{z) = A^i'i-A^z) + A 2 Bi'{-A 1/3 z). (120) 
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We must find solutions of this form which resemble the behaviour of the eigenfunctions of 
the equation with 7 > which obey the boundary conditions 

d2 { +(0) = and /_(0) = 0. (121) 
dz z 

The functions Ai' and Bi' have the following forms in the neighbourhood of z = 



Ai'(y) = ci(y + O(y )J -o 2 {l + Y + °(y )) ( 122 ) 

Bi'(y) = ^ Cl (^ + 0(y 5 ))+VSc 2 (l + ^ + 0(y 6 )) (123) 

with y = -A 1 / 3 ^, ci = 3" 2 / 3 /r(2/3), and c 2 = 3 -1 / 3 /r(l/3). So the positive-parity 
solution corresponds to the choice 

,4+ = _V3 and A+ = 1 (124) 

while 

AT/ = \/3 and = 1 . (125) 
At large values of y the corresponding wave functions are of the form 

f±{y)~{-y) 1/i ^( 1 ^-y? l2 + \±l)- i™) 

Noting that, near z = 0, f*dzp(z,\) = (2/3)(-y) 3 / 2 and (z/p(z, A)) 1 / 2 oc (-y) 1/4 , the 
solution coming from the origin is 

/o,±(*) = C7o,± [*M*, A)] 1/2 sin ( £dzp(z, A) + ~ ± |) . (127) 

The forms (j!17|) and l|127j) should be smoothly connected for < z < Zt. p .- This requires 

^ 2t ' P dzp(z,A) + |±^ = (n + l)vr, (128) 

for n = 0, 1, . . . together with C t . p . = (— l) n Co,±- Writing 

5(A) = r P 'dzp(z, A) = ^-(-4A + 2) , (129) 
Jo 11 



we find the quantisation condition 



S{\±) = (130) 



2 3 

Using (|129|) the quantisation condition takes the form 

A = / ~ 3n 6Ven P arit y (131) 
\ — 3ra — 2 odd parity 

which corresponds exactly to the spectrum (|2*9*)l obtained by exactly diagonalising H. 
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8 Results for other values of ( 



Up to now we have only considered the case of generic random forcing (where the constant 
D\ in equation (jHJ) is not zero), corresponding to the case £ = 1 in equation ®. In this 
section we explain two cases where other values of £ arise and briefly describe results for 
arbitrary positive values of C, analogous to the results obtained in section 6. 

First consider the case where the force is the gradient of a potential, f(x, t) = dV(x, t) jdx. 
We assume that the potential has mean value zero and correlation function C(X,T) = 
(V(x + X, t + T)V(x,t)} . Assuming that C(X,T) is sufficiently differentiable at T = 0, 
the diffusion constant is 



1 f°° I dV dV 
D(p) = - 2 j_Jt{^W m ,t)- M 

f°° d 2 C 
7 dX^(X,mX/ri 



2|p| 
— m 



oo 
oo 

dX 

-oo 



d 2 C , s m d 3 C , 



m 2 „„ d A C 



+ ^ X2 WXVT (X ' 0) + 0(X ' 



(132) 



Integration by parts shows that the integral over the first term of the expansion is zero, 
and the integral over the second term is zero by symmetry. The leading-order contribution 
in comes from the third term. Integrating this term by parts twice gives 

-to 3 f°° d 2 C 

m ~wL AX w^ o) - < i33) 

Thus in the case of a potential force with a sufficiently smooth correlation function we 
have C = 3. 

An exceptional case which is worthy of comment is when the potential V(x,t) is it- 
self generated from a set of Ornstein-Uhlenbeck processes Aj(t) by writing V(x,t) = 

Aj(t)<&j(x), where the &j(x) are elements of some suitable set of basis functions. In 
this case the correlation function of V(x, t) is of the form c{x) exp(— 7|t|) [for some function 
c(x)]. Then the second term in the expansion on the final line of equation (|132|) does not 
vanish by symmetry and we find D(p) oc \p\~ 2 , that is C = 2. 

For general positive values of £ the Hamiltonian ()22|) is replaced by 

2 4 1 1 dz\z\<dz v ' 

Its ground state 

A+ = and Vo + (^)=C + e- |2|C+2/(4+2C) (135) 

and first excited state 

Aq = -1 - C and = C o z \ z \ C e~W C+2 /( 4+2 ^ (136) 

are found by inspection. Raising and lowering operators can be introduced in a manner 
analogous to eqs. ([24126)1 . We write 

H = a~\z\-^a + (137) 
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with = d z ± z\z\*> /2. The operators 

A = a + \z\^d + and A + = a~|z|~ c a~ (138) 

satisfy 

[H, A] = (2 + C)i and [ff, i + ] = -(2 + C)i + . (139) 
and act as lowering and raising operators. For the spectrum of H we obtain 

A+ = -(2 + C)n and A" = -(2 + Qn - 1 - C . (140) 

These expressions replace (|29|) . Note also that the commutator of A and A + is 

[A,A+]=-(2 + Q(H + G) (141) 

where G = a + \z\~^a~ and H — G = I. The normalisation of the eigenstates 

i+IVO = C" +1 |C +1 ) (142) 
A^n) = C- \i>~_ x ) (143) 

is determined as in section 14.1.21 We obtain 

(C~ +1 ) 2 = (2 + C)(n + 1)[(2 + C)n + 3 + 2fl (144) 

and 

(C+ +1 ) 2 = (2 + C)(n + 1)[(2 + Qn + 1] . (145) 

The results of section [S] for the matrix elements Z mn = (^~|5|^ n ) and for ^(0)/^q (0) 
generalise as follows: 



+ r(^r -m + n) A /r(n + l)r(2^T + m) 

(-ir- m ^f— (m + n+1)— V (146) 

r (2+c) r(2 -m + n) Jr(|ty + n )r(m + l) 



^+(0)M+(0) - ( I)", / r[((2 + C)n + l)/(2 + Q] (147) 

vuuj/VoW - i y r(n+i)r(i/(2 + c)) ' 1 j 

This allows us to obtain, for example, the diffusion constant 



i 



1 (#B> ( \* (2 + Q-^fi B (^, J fe,^ ! ^,^ i i) 
"' 2 U 4+< / sm(^)r(|tf) 2 1 ' 

describing the dynamics at large times. Upon substituting £ = 0, eq. I)148J) reproduces 
the standard Ornstein-Uhlenbeck result, and it gives (|9*8*|) for £ = 1. 
For the short-time anomalous diffusion we obtain 

(x 2 (t))=C,((p^ c )W m - 2 ) i W (149) 
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with 



" ! /'■' . a(x) — a(xy) 



C x = -C dx x 2 +c / dy 
Jo Jo L 

where a(x) is the same as in (|113j) . and 



xo'(x) 



1 1 " (150) 



1+C , N ±K 

y 2+C (l _y)2 +C 



2(2 -ft) 2+c 

c = 1 y . 151 

This reproduces (|112|) for £ = 1 and concludes our summary of results for other than 
generic random forcing. 
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Figure 1: The spectrum of H consists of two equally spaced (ladder) spectra A~ and 
which are 'staggered' (that is, they are interleaved with un-even spacings). A and A + do 
not change the parity of the eigenfunctions. 
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z 

Figure 2: The trajectories of the classical Hamiltonian (|114l) are figure-of-eight orbits. 
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